Objectives: numerically simulate basic Markov chains (discrete time and discrete state space).
Setup: after retrieving the resources for the lab on moodle:
lab1 conda environment from the provided requirement.txt fileconda create --name=lab1 --file=requirement.txt
conda activate lab1
# do not forget to deactivate the environment if needed
# you can remove the environment once you are done
conda env remove --name=lab1
lab1)Assessment → grade /20 (possibly converted later on to a grade ranging from F to A (A+))
This lab session will be evaluated, based on your answer to the exercises reported in a Jupyter notebook (e.g., this one) and any additional .py file produced. In particular:
Additional evaluation criteria:
import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import package as pack
import pandas as pd
Consider a system of $K = 30$ particles (labeled from 1 to $K$) evolving in a closed box. The box is divided into two compartments in contact with each other, respectively identified by an index, $0$ and $1$. A hole at the interface between the two compartments allows the particles to move from one compartment to the other.
Particle motion is modeled as follows: at each discrete time intant $n$, one particle is chosen uniformly at random and moved from its current compartment to the other. Let $X(n)$ denote the number of particles in compartment $0$ at time $n$.
1. Briefly justify that $\bigl(X(n) \bigr)_{n \in \mathbb{N}}$ is a Markov chain.
Pour décrire une chaine de Markov, on décrit son squelette et les transitions entre les différents états.

Avec $q = 0$
Ainsi, le modèle est une chaine de Markov.
2. Is the chain irreducible? (Positive) recurrent? Aperiodic or periodic? Prove each statement from your answer (recall the main steps covered during the exercise session).
La chaîne est irréductible car $\forall i,j$ $ \exists chemin $ tel que $P(i->j)>0 $
La chaine est récurente car on a pour $T_0$, temps de premier retour: $$ P(T_0<+\infty) = \sum_{n=1}^{\infty}P(T_0=n) = 1 $$
La chaine est périodique de période 2 car la chaîne peut être décomposé en de sous partie sur son univers : nombre impaire et nombre pair et c'est le plus petit nombre pour décomposer l'espace des états comme ceci. Si $X_n$ est pair (resp impair), alors $X_{n+1}$ est impaire (resp pair).
3. Recall the structure of the transition matrix, and encode it in Python (without any for loop).
Hint: use the function
numpy.diag.
K = 30
p=1/2
diagonal = np.array([i/K for i in range(1,K+1)])
P = np.diag(diagonal,k=-1)+np.diag(1+1/K-diagonal,k=1) # Définit la matrice de transition
4. Numerically verify that the binomial distribution $\mathcal{B} (K, 1/2)$ is invariant for the chain. This invariant distribution will be denoted by $\pi$ in the rest of the exercise.
Hint: you can use
scipy.stats.binom.
B = binom(K, p)
x = range(K+1)
mu_bin = B.pmf(x)
plt.bar(x,mu_bin)
plt.title("Histogramme de la distribution binomiale")
Text(0.5, 1.0, 'Histogramme de la distribution binomiale')
On vérifie que µ est bien une mesure invariante avec la propriété : la mesure invariante est un vecteur propre à gauche pour P associé à la valeur propre 1
mu_bin.T@P - mu_bin
array([-5.16987883e-25, 3.30872245e-24, 7.94093388e-22, -3.81164826e-21,
3.38813179e-20, -1.08420217e-19, 1.08420217e-19, 6.50521303e-19,
-3.46944695e-18, -1.21430643e-17, 2.08166817e-17, 5.55111512e-17,
-1.38777878e-16, 1.94289029e-16, -2.22044605e-16, 3.05311332e-16,
-2.49800181e-16, 2.35922393e-16, -1.52655666e-16, 1.52655666e-16,
-1.04083409e-16, 3.64291930e-17, -4.33680869e-18, 1.73472348e-18,
0.00000000e+00, -2.71050543e-20, 3.72694497e-20, -1.27054942e-21,
1.16467030e-21, 4.30133919e-23, 2.27474668e-24])
On vérifie ainsi de $\mu_{bin}$ est une mesure invariante.
5. Implement a Python function ehrenfest to simulate a trajectory of the chain for a system of $K$ particles, for initial distribution $\mu$ (for instance a Dirac centered in $0$, meaning that the compartment $0$ is empty at $n = 0$). The maximum number of time steps will be controlled by an input parameter $n_{\max}$. Random number generation will be controlled by a random number generator passed as an input to the function.
For an efficient implementation, do not use vector-matrix product with the transition matrix implemented in 3: rely on the description of the system instead.
Cette fonction est implémentée dans le fichier package.py
6. Simulate a trajectory of the chain starting in state $0$ for $n_{\max} = 5000$. Display the result in function of the time index $n$. Briefly describe the curve you obtained.
simulate_ehrenfest = pack.ehrenfest(mu=0,K=30,Generator = np.random.default_rng(),nmax = 5000)
plt.plot(simulate_ehrenfest)
plt.title("Trajectoire de la chaîne")
plt.xlabel("index de temps")
plt.ylabel("Indice de l'état courant")
Text(0, 0.5, "Indice de l'état courant")
7. Compare the empirical histogram of the trajectory obtained in 5. to the theoretical limit distribution $\pi$. What do you observe?
hist = plt.hist(simulate_ehrenfest,density=True,bins=K-5)
plt.xlabel("Etats visités au cours de la simulation")
plt.ylabel("Probabilité empirique")
plt.title("Histogramme des états visités au cours de la simulation des urnes d'Ehrenfest")
Text(0.5, 1.0, "Histogramme des états visités au cours de la simulation des urnes d'Ehrenfest")
On observe à la limite la distribution semble décrire la distribution de la loi binomiale.
plt.hist(simulate_ehrenfest,density=True,bins=K-5,color="red",label="Histogramme de densité des valeurs simulées")
plt.bar(x,mu_bin,color = "grey",label="Histogramme loi binomiale")
plt.legend()
plt.show()
8. a) Modify the function defined in 1. so that it returns the return time to state 0, defined as $T_{0,0} = \inf \bigl\{ n > 0, X(n) = 0 \mid X(0) = 0 \bigr\}$.
print(f"Le temps de retour pour un K=10 est de (simulation sur un trajet): {pack.return_time_0(K=10,Generator = np.random.default_rng())}")
Le temps de retour pour un K=10 est de (simulation sur un trajet): 110
8. b) [Optional] Run several chains (about 5, ideally in parallel) for $K = 10$, $n_{\max} = 5000$, and compare the empirical average of $T_{0,0}$ to $\pi(0)$. What do you observe?
Hint: a good tutorial showing how to run functions in parallel in Python is available here.
T_list = [pack.return_time_0(K=10) for i in range(2000)]
plt.hist(T_list,bins=100,density=True)
plt.title("Temps de retour en 0 pour K=10 avec 2000 simulations")
plt.xlabel("Temps de retour en 0")
plt.ylabel("Probabilité de ce temps de retour")
plt.show()
L_k=[]
for K in range(1,20):
T_list = [pack.return_time_0(K) for i in range(2000)]
L_k.append(np.mean(T_list))
plt.plot(range(1,20),L_k)
plt.xlabel("Evolution du paramètre K")
plt.ylabel("temps de retour en 0 moyen ")
Text(0, 0.5, 'temps de retour en 0 moyen ')
8. c) Comment on the possibility of numerically observing the chain returning to its initial state as $K$ increases.
...
Let $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ be a discrete time homogeneous Markov chain defined by the following initial distribution $\mu$ and transition matrix $P$
$$ \mu = [0, 1, 0, 0, 0, 0], % \quad % P = \begin{pmatrix} 1/2 & 1/2 &0 &0 &0 &0 \\ 1/4 &0 &0 &1/4 &1/4 &1/4 \\ 1/2 &0 &0 &0 &0 &1/2 \\ 0 &1/2 &0 &0 &1/2 &0 \\ 0 &1/3 &0 &1/3 &0 &1/3 \\ 0 &1/3 &1/3 & 0 &1/3 &0 \end{pmatrix}. $$1. What can you say about the Markov chain $X$? (irreducibility, positive recurrence, periodicity, ...). Justify each of your claim, citing the relevant results from the lecture.
Ici les propriétés de la chaine de markov sont directes, nous détailleront peu les preuves :
La chaine est irréductible car il existe une unique classe d'équivalence notée $C_{eq} = {0,1,2}$ En effet pour chaque couple i,j d'états, $\exists$ un chemin de $i->j$ avec $P(i->j)>0$ Ainsi la probabilité d'aller de i en j en un nombre d'état fini est non nulle.
Cette classe d'équivalence est fermée donc elle est récurrente (une classe d'équivalence à soit tous ces états récurrents, soit tous ces états transitoires), donc tous les états sont récurrents.
Cette chaine est de période 1 car il n'existe pas de séparation de l'espace d'état en sous espace tel que la probabilité de passer d'un sous-espace à un autre soit de 1. (ie : Que l'on soit dans l'état 0, 1 ou 2 il est possible d'aller directement dans les deux autres états.)
2. Write a function simulate_dthmc simulating the trajectory of the Markov chain $\bigl( X(n) \bigr)_{n \in \mathbb{N}}$ for $n_{\max}$ time steps. The signature of the function should include the following elements:
To this end, do not use matrix vector products, which would lead to an extremely inefficient algorithm in this case.
P = np.array([
[0.2,0.7,0.1],
[0.9,0,0.1],
[0.2,0.8,0]])
distribution_init = [0,1,0]
3. Simulate a trajectory of the chain for $n_{\max} = 2000$ starting from $X(0) = 1$. Plot the histogram of the states visited by the chain.
nMax=2000
states = pack.simulate_dthmc(transitionMatrix=P,distribution=distribution_init,nMax=nMax)
states=pd.DataFrame(states.astype(int),columns=['Etat'])
pack.show_histrogramm(states=states)
4. Determine numerically an invariant distribution $\boldsymbol{\pi}$ of the chain (e.g., based on an eigendecomposition numpy.linalg.eig). Is it unique? Compare it to the histogram obtained in 2 (superimpose graphs). What can you conclude?
La chaine étant récurrente et l'espace d'état fini, nécessairement la mesure invariante est unique. En effet d'après le Th de Perron-Fröbenius : Pour $P$ matrice stochastique et irréductible, 1 est valeur propre de cette matrice $\exists !$ vecteur à coefficients positifs associé à cette valeur propre tel que $\sum_{i=1}^{S}\mu_{i} = 1$
w,v=np.linalg.eig(P.T)
pack.show_superImposed_Histogramms(states,v)
Vecteur propre associé à la valeur propre 1 : [-0.75529182 -0.64035611 -0.13956479] Vecteur propre associé à la valeur propre 1 normalisé : [0.49197861 0.4171123 0.09090909]
Nous constatons en effet que la probabilité de visite de l'état 2 est autour de 0.09 contrairement à celle de visiter l'état 1 ou 0 qui vaut 0.41 ou 0.49 Ainsi en simulant la chaine un grand nombre de fois, nous devrions tomber sur cette probabilité invariante.
5. a) Compute $\mu_n \triangleq \mu P^n$, the probability distribution of $X(n)$. What is the limit of $\mu_n$ as $n$ goes to $+\infty$? Illustrate the result numerically.
nStep=50
distribution_n,history = pack.compute_distrib_n(P,distribution_init,nstep=nStep)
pack.show_histogramm_vector_df(pack.convert_vector_to_df(distribution_n))
$\mu_{n}$ tend vers la distribution invariante de la chaine de Markov, l'histrogramme ci-dessus illustre bien cette convergence vers la probabilité invariante de la chaine.
5. b) Display on the same graph the curves $n \mapsto \mu_n(i)$ for $i = 1, \dotsc , 6$, and compare with $\pi$. Display on another graph the function $n \mapsto \Vert \mu_n - \pi \Vert_1$, where $\Vert \cdot \Vert_1$ is the $\ell_1$ norm. What does each of these curves illustrate?
history = np.array(history)
history= history.reshape(50,3)
plt.figure()
plt.plot(history[:,0],label="Evolution de la distribution en 0")
plt.plot(history[:,1],label = "Evolution de la distribution en 1")
plt.plot(history[:,2],label = "Evolution de la distribution en 2")
plt.xlabel("Nombre de steps")
plt.legend()
<matplotlib.legend.Legend at 0x15bd7df7ac0>
pi_vect = np.array([pack.norm_eigen_vector(v) for i in range(0,nStep)])
diff = np.linalg.norm(history - pi_vect,ord=1,axis=1)
plt.figure()
plt.plot(diff)
plt.title("Evolution de l'écart entre la mesure invariante et son estimation")
plt.xlabel("Nombre de steps")
plt.ylabel("Norme l1 de la différence")
plt.ylim(0,1)
plt.show()
Ces résultats montrent que la distribition des probabilités de chaque état converge ici très rapidement. (1er graphe) Le second graphe montre l'écart entre cette distribution théorique et celle estimée.
6. For each state $i \in \{1, \dotsc, 5 \}$, simulate 100 trajectories starting from the state $i$ until the return time to $i$. For each state, compute the (empirical) average return time. Compare with its theoretical value.
MoyennePerState=[]
for state_considered in range(0,3):
states = pack.simulate_dthmc(transitionMatrix=P,distribution=distribution_init,nMax=nMax)
positionOfReturn = []
for step in range(len(states)):
if states[step]==state_considered:
positionOfReturn.append(step)
MoyennePerState.append(pack.computeMoyenneReturn(positionOfReturn))
print(f"Le nombre moyen de retour pour chaque état est défini par (moyenne empirique)")
print(f"Moyenne du nombre de step pour le retour en 0 : {round(MoyennePerState[0],2)}")
print(f"Moyenne du nombre de step pour le retour en 1 : {round(MoyennePerState[1],2)}")
print(f"Moyenne du nombre de step pour le retour en 2 : {round(MoyennePerState[2],2)}")
Le nombre moyen de retour pour chaque état est défini par (moyenne empirique) Moyenne du nombre de step pour le retour en 0 : 2.04 Moyenne du nombre de step pour le retour en 1 : 2.36 Moyenne du nombre de step pour le retour en 2 : 11.27
En théorie : Pour une chaine de markov positive récurrente irréductible et apériodique $lim_{n \to +\infty}P(X_{n}=i) = \mu _{i} = \frac{1}{E(T_{i})}$ Avec µ la mesure invariante pour le système
On en déduit les valeurs théoriques des temps de retour :
tabTempsRetour = [1/val for val in pack.norm_eigen_vector(v)]
print(f"Moyenne du temps de retour en 0 : {round(tabTempsRetour[0],2)}")
print(f"Moyenne du temps de retour en 1 : {round(tabTempsRetour[1],2)}")
print(f"Moyenne du temps de retour en 2 : {round(tabTempsRetour[2],2)}")
Moyenne du temps de retour en 0 : 2.03 Moyenne du temps de retour en 1 : 2.4 Moyenne du temps de retour en 2 : 11.0
Il y a évidemment un écart du à notre estimation en temps fini mais il est clair que notre estimation converge bien vers le nombre théorique de retours.
import os
os.system('jupyter nbconvert --to html lab1_notebook.ipynb')
0